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Abstract 

A theory for multiple scattering of elastic waves is presented in a random 
medium bounded by two ideal free surfaces, whose horizontal size is infinite 
and whose transverse size is smaller than the mean free path of the waves. 
This geometry is relevant for seismic wave propagation in the Earth crust. 
We derive a time-dependent, quasi-2D radiative transfer equation, that de- 
scribes the coupling of the eigenmodes of the layer (surface Rayleigh waves, 
SH waves, and Lamb waves). Expressions are found that relate the small-scale 
fluctuations to the life time of the modes and to their coupling rates. We dis- 
cuss a diffusion approximation that simplifies the mathematics of this model 
significantly and which should apply at large lapse times. Finally coherent 
backscattering is studied within the quasi-2D radiative transfer equation for 
different source and detection configurations. 
46.40.Cd, 91.30.Dk, 43.20.Fn, 62.30.+d 
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I. INTRODUCTION 



Multiple scattering studies of elastic waves have become relevant to get to a deeper un- 
derstanding of the seismic Coda Q |T|fJ] and the principle of equipartition [JM. Multiple 
scattering is believed to occur in the spectral regime of 1- 15 Hz. Waves with larger fre- 
quencies suffer from absorption. Waves with smaller frequencies don't see the disorder any 
more. 

The wavelength of a transverse elastic wave at a frequency of 2 Hz is 1.7 km. This 
brings us to a basic problem in seismic studies of multiple scattering of elastic waves: all 
measurements take place at the free surface. As a result, they suffer from coherent reflections 
and mode-conversions. To be of any relevance to seismology, a multiple scattering theory 
should be capable of describing the boundaries in a very precise way, preferentially on the 
level of the wave equation. The standard equation of radiative transfer [Bj does not have 
this property. All phase information has been lost, and the accuracy of its spatial resolution 
is estimated to be some small fraction of the mean free path, estimated to be equal to 20 
km at least for seismic waves which is usually still much bigger than the wavelength. The 
equation of radiative transfer has been studied in ultrasonics for elastic waves at much higher 
frequencies , for which the complication of near- field detection is much less a problem. 

In this paper we present a multiple-scattering model that has been adapted to the needs 
of seismology. It incorporates the complex polarization properties of elastic waves as well 
as the mode conversions at the free surface. At the same time, we can formulate an almost 
convential radiative transfer equation describing mode conversions induced by scattering. 
Many other contemporary phenomena, such as equipartition and coherent backscattering, 
can be studied for measurements taking place near the free surface. 

The set-up of this paper is as follows. In section 2 we look at the wave equation for elastic 
waves, and we will define the Green's function for elastic wave propagation. In section 3 we 
introduce small-scale fluctuations and define the ensemble-averaged Green's function. This 
provides us with the extinction times of all elastic modes. They will serve to define our 
quasi-2D approximation. In section 4, the transport equation is derived, which describes 
the time-evolution of the ensemble-averaged energy-contents of all individual modes, and 
whose stationary solution exhibits equipartition of energy between all modes. In section 
5, we discuss the application of the well-known diffusion approximation to this quasi-2D 
model, introducing a N xN diffusion tensor for N modes. Finally, in section 6 we investigate 
coherent backscattering using our Quasi-2D approximation for different source and detection 
configurations. Section 7 is devoted to conclusions and perspectives. 
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II. TIME-EVOLUTION OF ELASTIC WAVES 



In this section we will formulate the mathematics of elastic wave propagation in a way 
that is suited to apply conventional methods in wave transport. Many elements have already 
been discussed very thoroughly by Papanicolaou et al. 0, and some will be recalled here for 
convenience. We start out with Newton's second law for the elastic displacement u at time 
t and position r, 



p{r)d 2 t Ui = djGi^r) + fi{r,t) . 



(1) 



Here, p(r) is the local mass density, f(r,t) is an external force per unit volume, and a ij(r) 
is the local traction which, by Hooke's law, is given by [p^0| -[T2| 



Cijki(r)£ki{ T ) 
\{r)E kk 5 i:j + 2/i(r)£jj(r) 



(2) 



with Ski = \{dkUi + diUk) the stress tensor. As always, summation over repeated indices is 
assumed implicitly. The second equality applies to an isotropic elastic medium, in which 
case the four-rank tensor Cyjy can only have two independent contributions, proportional 
to the Lame moduli A and fi. Following Papanicolaou et al. 0] we shall introduce the vector 
field, 



I slh ■ u \ 



I 

Ud t Ui 



(3) 



This vector has 9 independent components since £jj is a symmetric tensor, whose trace is 
given by the first component. We have introduced the operator p = — iV. It can readily be 
checked that Eqs. (HD and (0) combine to the following time-evolution problem, 



i0 t |*(r,t)) = K(r,p) -|*(t)> + |*/(t)> 



(4) 



with the time-evolution operator, 

/ 



K 



V on 



\ 

5 II / 



(5) 



and the external force-term ^j(r,t) = (0, — f (r, t)/yp(r), 0). We have introduced the third 
rank tensor /^(p) = | (pi8jk + p5ik) and used the formal Dirac notation for vector fields 
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to facilitate later the more convenient mode base. The number of arrows determines the 
order of the tensor. For clarity, we have put horizontal arrows when they contract in a right- 
hand side product with a vector. If A and /i are real-valued, the matrix K is manifestly 
symmetrical with respect to the ordinary Cartesian scalar product. As a result, 

(*(t) | *(t)) = J dr*(r,t)* ■ *(r,t) = J dr^A(r)(V • u) 2 + ^p(r)(«9 t u) 2 + /i(r)Tr e* ■ e, 

(6) 



recognized as the total elastic energy [|TJ], is conserved in time if no external forces are 
present. It is customary to split off a term \x (curlu) 2 (describing pure shear wave energy) 
and 2/x (div u) 2 (contributing to compressional energy) from the last term, leaving a rest term 
I. This identifies four terms as "kinetic energy", "compressional energy", "shear energy" 
and an interference term |T^] . The latter vanishes for plane waves with either pure transverse 
or pure longitudinal polarization. 

Equation (f|) can easily be Laplace-transformed (Im z > 0). This yields the solution 

|*(z)) = [ 2( -K]- 1 [i|*(t = 0)> + |*/^)>] • (7) 

The operator [z — K] _1 = G(z) will be called the Green's function, and is introduced here 
for future need. It is convenient to define t — just before the source sets in so that 
vE^i = 0) =0 and the force field becomes the source for wave propagation. 



III. EXTINCTION OF ELASTIC WAVES IN A LAYER 



We consider a homogeneous elastic plate with mass density p of infinite horizontal di- 
mension and of thickness H. Both sides of the plate will here be assumed to be free surfaces, 
with traction-free boundary conditions. The Lame coefficients are A and fi, in terms of which 
the P-wave speed is a = \J (2/x + A)/p and the S'-wave speed is (3 = \ f\ifp- The eigenmodes 
have been discussed in great detail by Weaver [|H|[14]]. Their representation (Q) can be ob- 
tained straightforwardly and we shall denote them by {^n}. The index n is a discrete index 
that labels, at constant frequency, the symmetric and anti-symmetric Lamb and SH waves 
in the plate, including the two Rayleigh surface waves (Fig. |I]). 

The structure of the symmetric and antisymmetric Lamb modes is extremely rich whereas 
the SH modes are "simple" guided waves. For the sake of clarity let us focus on the 
antisymmetric branches, indicated with normal lines in Figure |I[ The first antisymmetric 
mode (first black dot on the right in Figure [TJ) is an antisymmetric Rayleigh surface mode. Its 
displacement is evanescent for both the compressional and the shear component. Rayleigh 



4 



modes propagate somewhat slower than bulk S or P waves. As a result they lie on the 
right of the two dashed lines indicating the pure shear and pure compressional excitations. 
The second antisymmetric Lamb mode (third black dot) lies between the two dashed lines 
indicating the pure shear and pure compressional excitations. This mode is evanescent for 
its compressional component but has a propagating shear displacement. It behaves like a 
pure shear mode as we go away from either one of the free surfaces. As a matter of fact 
its potential energy is mostly shear since its compressional potential energy is negligible. 
Finally the antisymmetric mode most on the left in Figure [3] lies on the left of the lines that 
indicates pure shear and pure compressional excitation. Even deep in the plate this mode is 
a mixture of P and S displacements. As we increase the frequency, the organization of Lamb 
modes stays intact. One finds two surface Rayleigh modes (symmetric and antisymmetric), 
"evanescent P but bulk S" modes and modes that are both bulk S and bulk P. 

By translational symmetry, the eigenmodes can be chosen proportional to transverse 
plane waves with wave number k. We will treat him initially as exp(zk ■ x)/y/A, with a 
discrete contribution of k to the label n as a result of the periodic boundary conditions on 
both sides of a square plate with surface A, and eventually take the limit A — > oo. 

We will now assume the presence of disorder in the plate, to be specified more precisely 
later on, and calculate the Green's function, averaged over this random disorder. The exact 
meaning of this averaging in seismic observations will be addressed elsewhere. This procedure 
is the first part to formulate a transport theory [|Tj|. Let the disorder be represented by 
a perturbation <5K in the time-evolution operator: K = Ko + <5K. The ensemble-averaged 
"retarded" (outgoing) Green's function at frequency u is given by, 

This "Dyson" equation defines the mass-operator X(a>). The lowest order contribution is 
given by |T5| ], 

S(cu) = /5K 1 ■ 5k) + 0{5Kf . (9) 



\ u + iO - K / 
Next, we can insert the complete and orthonormal set {*f? n } of the homogeneous plate, 
defined above. Standard first-order perturbation theory yields, 



with 



SnM = £ (|(*„| 5K |* m )| 2 ) L— (11) 
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The imaginary part of this parameter is negative, and is identified with — l/2r n , where r n is 
the extinction time of mode n. 

In general, both p(r), A(r) and /x(r) are random variables. We will simplify the problem 
by assuming that p(r) is constant, and that velocity fluctuations are due to fluctuations in 
the two Lame coefficients: A(r) = A + SX(r) and /x(r) = fi + <5/i(r), with A and fj, the 
coefficients of the homogeneous layer. In that case, 



5K= — 

VP 



(5A(r)/2v^)p 
01 



p | (5X(v)/2VXo) 

0|| (Mr)/2v / 2^)L(p)|| 





£(p) | (^(r)/2v^) 

8u 



(12) 



A straightforward calculation, employing integration by parts, finally leads to, 



\(*n\5K\% 



UJ 



dr J dr' 

{ (a(r)5A(r')> (V • u n )*(V • u m )(V • <)*(V • 
+ (5Mr)5Mr'))Tr<-, m Tr(4r-4 
+ (5A(r)<^(r')> (V • u n )*(V • u m )Tr (e'J* • 4 + cc. } 



(13) 



To evaluate H n (u) we must specify the spatial correlations between the Lame coefficients. 
The simplest choice is to assume correlations that are short range with respect to the wave- 
length, 

(5X(r)SX(r')) = a 2 x (z)S(r - r') (14a) 

= - r ') (14b) 

(5pL(r)5X(r>)) = ^ A (^(r - r') . (14c) 

Without extra difficulty, we can still allow for a depth dependence of the correlation func- 
tions. S n can now be evaluated for a big plate for which J2m -^HiA J d 2 k/ (27r) 2 , including 
a sum over the different branches. All factors A cancel if a transverse plane wave normal- 
ization exp(ik • x) is adopted. For the extinction time of mode branch j at frequency uj, we 
find 



1 



UJ' 



d 2 kj 

2tt 



(15) 



with riiiyj) = ki{uj) / Vi{uj) in terms of the group velocity Vj = dcUj/dk;. The "mode scattering 
cross-section" is defined as, 

2 



|V • Uik; I |V • Uj kj | 



+ 2al x (z)Re (v • u* k; V • u jkj Tr £ * ki • e jkj ) } 



(16) 
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We have chosen to split of the factor rij, so that this matrix is symmetric. According to 
our model the extinction time Tj does not depend on the direction of the horizontal wave 



number kj. 



The imaginary part of the ensemble-averaged Green's function is directly related to the 
excitations of the waves [|17]]. The spectral density J\f(uj) per unit surface can be expressed 



as, 

A^) = -i-TrImGH = -Y: / T 1Z^__ . (17) 

V 7 ttA v ' tt^J (2tt) 2 [lu-lu^ + I/AtI v ' 

Due to scattering, all modes are spectrally broadened. The separation in wavenumber of 
two adjacent modes with the same frequency (see Figure [I]) is typically of order 1/H. The 
uncertainty in k is typically l/t>ik r ik, with the group velocity of the mode. If 

T ik >H/v ik , (18) 

one can assume that different modes at fixed k do not overlap, except at a few degener- 
ation points where the dispersion curves for modes with different symmetry (i.e. SH and 
Lamb) cross. This assumption is the Quasi Two- Dimensional Approximation (Q2DA). Cri- 
terion (|18D is typically satisfied in the Earth crust, which has H pa 30 km, a typical wave 
speed (3 ~ 3.5 km/s and a mean free time r > 15 sec. In the Q2DA we find for the spec- 
tral density per unit surface N(uj) = (27r) _1 J2i n ii showing that n i: defined in Eq. (|15D, 
represents the spectral weight per unit surface of mode i at frequency u> in phase space. 

In the following, all time scales will be normalized by the mean free time of S waves in 
an infinite medium with the same amount of disorder. This time depends only on a 2 which 
can be related to the correlation length and the shear velocity fluctuations. For a velocity 
fluctuation of 2% and a correlation length of 700m, both being typical seismic values, we 
get a shear mean free path pa 15s. Note that the correlation length is much smaller than 
the wavelength X s = 1700m, which justifies the short range approximation of Eqs (|14cj) . 

Figure |2] shows extinction times for different modes index, calculated from Eqs. (|15"D 
and fllB]), normalized by the mean free time of 5*- waves in an infinite medium. The plate 
thickness is H = 20.2As, which has iV = 106 modes. The disorder is chosen to be uniform 
in the whole plate, and the spatial correlations among the Lame coefficients is taken equal: 
a\ = a 2 = cr 2 A . SH modes show an extinction time very similar to the extinction time of 
S- waves in an infinite medium r s °°. On the other hand the Lamb modes present a more 
complex pattern: Rayleigh modes clearly show a shorter extinction time, Lamb modes with 
an evanescent compressional component behave very much like a bulk S wave. Finally, Lamb 
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modes with both bulk compressional and bulk shear components behave in a complicated 
fashion but tend to have an extinction time larger than S waves in an infinite medium. 

In the case of dominant ji correlation, a\ <C a 2 (dominant shear velocity fluctuations) 
the Lamb modes with both "bulk" compressional and shear components will have a some- 
what relatively larger extinction time. On the other hand, if the A correlation dominates, 
°a ^ a % (strong compressional velocity fluctuations), the same Lamb modes with 'bulk" 
compressional and shear displacements will have the shortest extinction time. 

We would like to point out that the life time of Rayleigh waves is not well described by 
our model since they suffer most from surface disorder (fluctuations in height), which is not 
included in Eqs. ([14]). 



IV. TRANSPORT EQUATION IN A LAYER 



The next task is the formulation of an elastic transport equation in the Quasi 2D 
Approximation. Basic observable is the ensemble-averaged intensity Green's function 
(G(u~) x G(u + )*), with u; ± = uj ± ^Q. It can be expressed in the complete base {*& n } 
of the homogeneous plate, giving rise to the matrix element L(u>, fl) nn ' mm ' (Figure 2). The 
Bethe-Salpeter equation [[I5|,[nj for this object reads, 



£-nn'mm'{^i^) ~ G n (uJ + )G n > (to )* 



SnmSrim' + X/ ^nn'll'i^, ^) £-11' mm' 'j ^) 
IV 



(19) 



with G n the Dyson Green's function defined in Eq. (§), and a new object U called the 
Irreducible Vertex. Upon introducing AG nn i(u>,fl) = G n (uj + ) —G* n ,{u)~) (idem for AS) this 
equation can be re-arranged into, 



[Q + {u n - a;*,) - AE nn /] C r 



(u,Q) — AG nn '(u,Q) 5 nm 5 n i m i + 

IV 



(20) 



This equation is still exact. We will now carry through a number of approximations relevant 
to our problem. For small disorder, the vertex U is given by [Tj| . 



(21) 



For short-range correlations, as specified in Eqs. (0), the vertex U can be straightforwardly 
related to the cross-section W(iki, jkj) defined in Eq. (p!|) . For typical wave packets is 
fl < w (i.e. a wave packet contains many cycles) so that we neglect Q in any functional 
dependence on frequency ( "slowly varying envelope approximation" ). The index n consists of 



S 



one discrete branch index j, and one index k that becomes continuous as A — > oo. The Q2DA 
neglects all overlaps between different branches, so that AG(u, £l) nn > — > 2mSjji6[u — Wj(k)]. 
If we let k — k' = q, and S m (uj) the source in mode representation, a new observable quantity 
Ljk can be defined as 

£nn'mm'(u, Q)S m S^, = 2tt5[uJ - LJj^Sjf x L jk (q, Q). (22) 

mm' 

In space-time the Q2D transport equation reads, 

L ik (x,i) = \S jk (uj)\ 2 8(t)5(x) 

W£ / ^n^O^-,/^)^ (x,t). (23) 

We will use this equation as a starting point for our calculations. The equation is essentially 
two-dimensional, with a finite number of modes (of order 2Huo/ (3) to take care of the third, 
vertical dimension. The great advantage of this equation is that the boundary conditions of 
the elastic waves have been dealt with exactly, i.e. on the level of the wave equation, contrary 
to conventional transport equations 0Hp|. We see that Lj k (x, t) can be interpreted as the 
specific intensity of the mode (j'kf) at frequency u, at horizontal position x, at a time t after 
the release of energy by the source. 
The source term Sjk(uj) is given by, 

Sfrbv) = (¥ jlt | 9 f )=uJ d 3 rf(r, W ) ■ u jk (r) . (24) 

Since Uj k is an eigenfunction for which the energy (^) has been normalized, we see that |Sjk| 2 
has the dimension of energy. Since (f2, q)-dependence has been neglected in the source, it 
emerges in our transport equation as a S(t)8(r) in space-time. 

A. Dynamics of the Equipartition Process 

Equation (p3|) has one very important property that has been discussed in great detail in 
the literature. By recalling the expression (|X3|) for the extinction time, it follows immediately 
that the specific intensity with the property that its total mode energy 

J d 2 xLjkj(x, t) = constant, (25) 

is independent of the mode-index j and independent of the horizontal direction of propaga- 
tion k, is a stationary solution for t > of the transport equation. All solutions converges 



d t + Vj ■ v + — 
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to this solution regardless the nature and position of the source. This implies that finally 
all modes have an equal share in the total energy contents of the plate. This phenomenon 
is called equipartition ^,W[-21\, and is believed to be a fundamental feature of the solution 
of most transport equations at large lapse times, provided absorption is absent, or at least 
small p^| . According to our definition (|22|), the total spectral energy per unit surface in the 
regime of equipartition is given by, 

d 2 k 



E »$) = E/ 7^2 / d ' x L jk (x,t)27r^-^ jk ) 



(2tt) 2 
constant x }^ rij . 



(26) 



We will introduce the spectral energy density Ei(x,t) of mode % per unit surface, and its 
current density Jj(x, t) according to, 



^ 2 k 

(2 

d 2 k 



— — - 2ir5 (u - Ufo) L t -k(x, t) = rti / — — L ik .(x, t) , 

[III) J Z7T 



/cl7k f d 2 k 

— — 2n5 (u - uJiv) ViL ik (x, t) = m J — v i Ir iki (x, t) 



(2 

An exact equation of continuity can be found from Eq. (^3[) by integrating over kj 

d 2 k 



(27a) 
(27b) 



dtEifc t) + V- Ji(x,t) 



2tt 



I'S'ikiH 



5(x)5(t)-^C t ^(x,t) 



(28) 



with the "mode-conversion matrix" . 



r - 6i i 

^ij — 

Ti 



d 2 ^ 

2tt 



(29) 



The mode-conversion matrix C has an eigenvalue with (left-hand) eigenvector {rij}, asso- 
ciated with the equipartition. The N — 1 nonzero eigenvalues, which can be called "Stokes 
parameters", of the mode conversion matrix C determine the dynamics of the equiparti- 
tion process. It depends on the initial conditions, i.e. how the initial release of energy is 
distributed among the different modes, as described by S^u). 

Figure |3] shows all eigenvalues of the matrix C in the case of a plate of thickness H = 
2O.2A5, for which the number of modes is N = 106. The disorder is uniform in the whole plate 
and the spatial correlation among all Lame coefficients is chosen equal: o\ = cr 2 = er 2 A . The 
time scale has been normalized to the mean free time of S-waves in an the infinite medium, 
which has the same amount of disorder, i.e. as described by Eqs. (|H[) . 

The largest eigenvalue (associated with the shortest life-time) has an eigenvector made 
of the symmetric and antisymmetric Rayleigh modes. This configuration is very sensitive 
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to the location of the disorder in the plate. If the plate does not have any disorder close 
to the two free surfaces (at the length scale of a wavelength) the Rayleigh modes, which 
have a penetration length of the order a wavelength, do not suffer from the disorder. As 
a consequence, their life-time would become very large compared to the mean free time 
of ^-waves in an infinite medium. On the other hand, if the disorder is localized close to 
the free surface the Rayleigh modes end up with a very large eigenvalue. The eigenvectors 
associated with the flat plateau in Figure [5] consist of modes whose shear component strongly 
dominates over the compressional part. As a result their eigenvalues are very similar to the 
inverse shear mean free time of a S wave in an infinite medium. The eigenvectors associated 
with the eigenvalues smaller than unity exhibit a strong compressional component. They 
are associated with longer life times as shown in Figure Q Quite logically they show up with 
a smaller eigenvalue (associated to a longer life-time) in the mode-conversion matrix C 

In the case of dominant \i correlation, o\ <C a 2 , the picture does not change drastically 
since Lamb modes are always dominated by shear. For dominant A correlation, o\ ^> 
a 2 the structure of eigenvalues of the mode-conversion matrix C is modified considerably. 
Eigenvalues that were previously associated with "bulk" P and S vectors now see their life- 
time becoming much shorter. An eigenvector with a small eigenvalue in Figure achieves 
a large eigenvalue. 

Figures [| show, for different kind of sources, how the initial release of energy is distributed 
among the different modes. Figure shows an isotropic explosion at a depth A s /3 from 
the free surface. An explosion is a purely compressional source, and does not excite any SH 
modes. Among the Lamb modes it excites preferentially the modes that are "bulk" for both 
compressional and shear components as well as Rayleigh modes. A source at a larger depth 
will not excite the Rayleigh modes since they have a penetration length of the order of the 
wavelength. 

Figure £|6 applies for a double couple in the xy plane at a depth A s /3 from the free surface. 
Contrary to the isotropic explosion, the double couple in the xy plane strongly excites the 
SH modes. Since the source is close to the free surface Rayleigh modes are excited as well. 
The Lamb modes which are "bulk" for the shear component but only evanescent for the 
compressional component are also excited. 

Figures f|c, d show the mode distribution for a double couple in the plane xz for two 
different depths of the source, A s /3 and 5A S . When the source is located close to the free 
surface the majority of the energy is distributed among the Rayleigh modes. Two Rayleigh 
modes are out of scale in Figure |]c and carry half of the released energy. Alone they carry 
half of the total energy released. On the other hand, when the source is situated deep in the 
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plate the pattern becomes very rich. One can see that the Rayleigh modes are not excited 
anymore. 



V. DIFFUSION APPROXIMATION 

Despite the many simplifications that have been carried out, the final transport equation 
(f23|) is still difficult to solve numerically. In future work, we intend to adapt our Monte-Carlo 
simulations, developed to solve the 3D radiative transfer equation fl2"3"| , |2"3]| , to this modified 
equation. In this section we shall carry out a final and rather familiar simplification, that 
facilitates a numerical solution. 

The diffusion approximation is typically valid at large lapse times, when currents start 
to become small. In that case, the specific intensity of mode i can be written as, 



Hi 



2 

E l (q,n) + — Vi - Ji(q,fl) + - 

vf 



(30) 



with rii = ki/vi the density of mode i in phase space introduced earlier. In real space q 
transforms into the 2D gradient V. Inserting the series (|30|) into Eq. (^) leads to the 
relation 

J i (r,0 = -EAiVE 7 -(r,t). (31) 



This relation is recognized as a generalized Fick's Law [25], generalized, because it involves 



different individual modes at the cost of one dimension. The diffusion matrix is given by, 

(D-').. = 2(^-^/^ H /( s k„ ! k j )^) , (32) 
v >*3 \vfTi rij J 2tc vfvj J 

It is easy to check the following relation, 

^ = — . (33) 
Dji rij 

Combining Eqs. (|31]) and (|28|) and transforming back to space-time yields the generalized 
2D diffusion equation, 

dtEfat) - J2 A,HV 2 ^(r,t) = Si(w)8(t) - E^H£ 3 (r,t). (34) 
i 3 

This diffusion equation is an ordinary partial differential equation that can be solved by 
conventional means. For an infinite plate no boundary conditions have to be specified: the 
boundary conditions at the two free surfaces have been taken care of exactly. For this reason, 
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the Q2D diffusion approximation is not expected to break down near the boundaries, as was 
noticed by Turner and Weaver for the conventional diffusion approximation 0. 

Equation fl3H]) still captures the time-evolution of the different elastic modes of the plate, 
and can thus be used to study polarization properties. Integrating equation (j34l ) over the 
horizontal coordinate r gives for the time evolution of the total modal energy 

dtEi(t) = Si(u)6(t) - J2 CiM E i$)- (35) 

3 

In fact this equation follows directly from Eq. (|28|) without the need to apply the diffusion 
approximation. Its formal solution is Ei(t) = J2j [ ex P ( — Ct)]^- Sj(u)6(t > 0). This can 
easily be evaluated using the complete set of eigenmodes of C. 

Figure |5|a shows the time evolution of the energy among the different modes for an 
isotropic explosion at a depth A s /3 from the free surface. The initial modal energy distribu- 
tion was shown in Figure [pa. For the sake a clarity we only display the evolution of three 
sub-classes of modes (Rayleigh, Lamb, SH) and not the whole distribution. Rayleigh modes 
are excited but not SH modes since the source is purely compressional. As time goes on, the 
mode occupation changes as a result of the dynamics of the equipartition process and finally 
tend to the equipartitioned distribution which does not depend on nature and location of 
the source. 

Figure |5|6 shows the time evolution of two "observable" energy ratios measured at the 
free surface: the ratio of shear to compressional potential energy, E s /E p , and the one of the 
horizontal to vertical kinetic energy H 2 /V 2 . After a few shear wave mean free times, the 
energy ratios stabilize to their predicted equipartition value E S /E P = 7.19, H 2 /V 2 = 1.77 
||. The ratios E s /Ep and H 2 /V 2 increase monotonically which is due to the compressional 
nature of the source. 

Figures |5]c, d present the equipartition process for a double couple source deep in the 
plate (5 A s from the free surface). For such a source the Rayleigh modes are not excited 
while the other Lamb modes and SH modes are strongly excited (see Figure |](i). The initial 
ratio of shear to compressional energy at the free surface is higher than the one for the 
explosion source due to the shear nature of the source. However, in both cases the energy 
distributions converge towards an equipartitioned distribution which is independent of the 
nature of the source and its location. Note that, for an exploding source, the equipartition 
process takes much longer a time, typically 6r^°. For the double-couple source in Figures ||c, 
d it is typically r s °°. 

It is not very difficult to show that in the equipartition regime, the generalized diffusion 
equation (|34]) further simplifies to a genuine 2D diffusion equation for the total energy 
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density, 

d t E(r, t) - D(uj)V 2 E(r, t) = S(u)5(r)6(t) , (36) 
with diffusion constant, 

Eq^M^M (3?) 

and source, 

%)=&/^ |S iki H| 2 . (38) 

Equation (|37| ) is recognized as an equipartitioned sum of all diffusion matrix elements. A 
similar result was obtained for the diffusion constant in an infinite elastic medium, in terms 
of the individual matrix elements for P and S waves [f7,21"],2l)f. Equation ([36]) has the simple 
solution, 

i.e. the local energy basically varies as E(u) ~ t^ 1 x S{uj) / D{uj) at large times. 

Table | shows the evolution of the ratio D(u) / D°°{uj) as a function of the number of 
modes in the plate. D°°(ijj) is the infinite medium elastic diffusion constant, obtained by 
Weaver [0] and Ryzhik [p7fl , with the same amount of disorder, i.e., as described by Eqs. ([14]) . 
The ratio changes form 0.72 for iV = 3 modes to 0.85 for a thick plate, i.e. H m I* . Our 
quasi-2D approximation starts to break down when the thickness of the plate exceeds the 
mean free path. 

VI. COHERENT BACKSCATTERING NEAR THE FREE SURFACE 

Coherent backscattering of waves is an interference effect that survives multiple scatter- 
ing. It refers to a coherent enhancement of intensity near the source ||26|| . The effect has 
recently been observed with acoustic PS| and elastic waves PU|PU . 



We recently investigated coherent backscattering of waves in a more seismic context 
32| , |3"3"|| . Specific aspects such as symmetry of the source, near field, leakage and the exact 



measurement process have to be understood before any seismic experiment can be consid- 
ered. Our analyses so far have been done either with scalar (acoustic) waves in a disordered 
plate with leakage or with fully elastic waves in an infinite medium [RBI . The last study 
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established that the enhancement factor of coherent backscattering is highly dependent on 
both the nature of the source and on the precise parameter that is being measured. More 
specifically, a measurement of simply (ui(u) 2 ) of waves released by a "double-couple" source 
will hardly give rise to a coherent enhancement, so that observation is unlikely. On the 
other hand, the measurement of (divu(a;) 2 ) of waves released by an explosion source maps 
exactly onto the acoustic problem, which has the maximal enhancement factor of 2. 

Both approaches are unable to model the coherent backscattering effect of wave propaga- 
tion in the crust, whose elastic eigenmodes are not plane waves. In addition, a measurement 
necessarily takes place at the Earth surface, whereas the source (an earthquake or explosion) 
can be located at depth. In this section we will investigate coherent backscattering using our 
Quasi 2D transport model. Recently, De Rosny et al. and Weaver et al. [^] reported 
the studies of coherent backscattering of elastic waves at frequencies around 1 MHz. 

Our analysis will closely follow the one given in Ref . . Starting point is the calculation 
of the vertex L nn / mm /(k, k', q) defined in Eq. ( p2|) and describing the ensemble-averaged, 
incoherent scattering of the modes (i, k+|q) and (i', k— |q) into (j, k' + ^q) and (j 1 , k' — ^q). 
By the reciprocity principle this object must be symmetrical with respect to left and right- 
hand indices. The diffusion approximation, applied to our Q2DA model yields for large lapse 
times, 

L Wjf (k, k , q) = — - - u/Q . (40) 

An inverse Fourier transform with respect to Q provides the time-dependence of the envelope 
of a wave packet with central frequency uo []. Similarly, the spatial dependence is obtained 
by an inverse Fourier transform over q, k and k'. The result is, 

L^y(u;,t,xi,X2 ^x 3 ,x 4 ) = ^ ^ ^ Su'Sjj' TiiUj J (kiX 12 ) Jq (kjX U ) ■ (41) 

The depth (i.e. z) dependence can by obtained by summing over the N eigenfunctions 
*&i(z) at frequency u. The coherent backscattering is due to constructive interference of 
time-reversed waves. It can be constructed straightforwardly by interchanging the indices 
(i'x 2 ) and (j'x 4 ) 



1 Seismic measurements usually have access to time-correlations (ipj(t — ^T")^^ + 2 r )/ °^ differ- 
ent components of the wave function (p3[). The smooth time-dependence of the envelope of the pulse 
at frequency u can be obtained by Fourier transforming this object with respect to r. Signal-to- 
noise can usually be increased by averaging over a time window At and using a finite bandwidth 
Auj over which the signal is not expected to vary too much. 
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Cuijfiuj, t, xi, x 2 -> x 3 , x 4 ) = 6XP ^ D ~T~^~ n i n i J o ( k i x u) Jo (kjX 32 ) . (42) 



Both L and C contribute to (G(z,xi — > j, x 2 )G*(i', x 3 — > j',x 4 )), but C survives only close 
to the source, as we shall see. To calculate actual enhancement factors, we must specify 
source and detector. In Eq. ( [24]) the source was already expressed in the eigenmodes (J, k). 
Different sources will now be considered. 



A. Monopolar source at depth 

We consider the source f(r) ~ io(uj)S^ 3 \r — ro), which represents a highly directional 
force field at position r , small compared to the wavelength. Equation (|4j) gives Sj^(u) ~ 
wfo(w) • Ujk.(^o) with zq the depth of the source. To simplify the analysis we will assume 
that the force is directed into the z-direction. This configuration was also studied by De 
Rosny et al. [^,0] using a thin chaotic 2D silicon cavity, with only 3 excited Lamb waves. 
In addition, their detection method of heterodyne laser interferometry is only sensitive to 
the normal displacement u z (z = 0). In seismology, the force field above may be a simple 
model for a volcano eruption. 

Let x be the horizontal distance between source and detector. The measured "incoherent" 
background is found from Eq. (^T|). 

l{x, t) ~ exp( ~f /g) /oM 2 E * k,*(°)l 2 E "i k^o)l 2 , (43) 

which is independent of x, but still depends on the depth zq of the source. The "coherent" 
contribution follows from Eq. (03) 



C{x, t) — f (u) 



E n * Ui, z (0) u itZ (z )* J (kix] 



(44) 



As was already mentioned in previous work, the ratio (L + C)/L, the so-called "enhance- 
ment factor", is independent of time at large lapse times J32J. An application of Cauchy's 
inequality shows that (L + C)/L < 2, with equality only if x = and if Ui tZ (0) = Ui tZ (zo) 
for all modes i. This can only be true if Zq = i.e. the source must be near the surface. In 
practice, to produce any measurable enhancement factor, the source must be at a depth less 
than the typical wavelength, as shown in Figures |B|a and ^p. A source with a force direction 
different from normal would have a lower enhancement as well. Note that the enhancement 
is symmetric in azimuth around the source. 
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B. Isotropic Explosion 



An isotropic explosion at depth zq is described by the force field f (r, uj) = B(u)V5(r — r ) 
|TT| . It can easily be shown that Sik(u) = — B{uj)uj divu i k(co | ). For a fixed frequency this 
depends on the mode label i but, very conveniently, not on the direction k of horizontal 
propagation. 

Let us first suppose that we measure the normal component of the displacement vector 
at the surface. Incoherent background and coherent enhancement are given by, 

exp(— ut/Q) 



L(x,t) 



Dt 



Hi Ui 



J2 n i l divu j( z o) 



C(M)~ ^f /Q) B{uf 



Ui M ii2 (0)div Ui(z )* Jo(kix) 



(45) 



The resulting enhancement factor (L + C)/L is plotted in dashed lines in Figure [?]a as 
a function of the horizontal distance, and in Figure ffi for a measurement on top of an 
explosion source as a function of the depth zq. Note that the enhancement never reaches 
its maximum value 2, not even when z = 0. In an infinite medium, a measurement of any 
component of the displacement vector of waves released by an explosion source would have 



had no enhancement at all near the source ||33|| . Here, the finite enhancement is due to the 
nearness of a free surface. 

The enhancement factor can be restored by a measurement of the dilatation (div u) in 
which case, 



L(x,t) 



exp(—ujt/Q) 



Dt 



B(uj) 2 ^2rii |divui(0)| 2 |divuj(z ) 



C(x,t) ~ B(uo 



2 exp(—ut/Q) 



Dt 



E riidw Ui(0)div Ui(z )* Jo(k;x) 



(46) 



A measurement of the dilatation restores the symmetry between detector and source, and 
reveals the maximum enhancement factor 2 when the detector is located close to the source 
as shown in solid lines in Figures [7|a, 6. 



C. Dipolar Source 

We next consider a single couple at the surface with normal displacement vector, and axis 
along the x-axis. This source can be represented by the dipole f (r, t) ~ d(uj) z^^^^r — r ). 
Such a source can be generated with laser interferometry on an elastic plate, and the resulting 



coherent backscattering effect was recently studied experimentally by De Rosny et al. pi ]. 
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The spatial derivative in the source finds its way in the Bessel functions, in the same way 



as was done in earlier work for the infinite system [ 3j| . We derive, again for a measurement 
of the displacement vector in the direction normal to the surface, 

exp(— tot/Q) 1 



L(x,i) 



Dt 



m Ui 



Y^rij \u jtZ {zo)\ k 



C(x, t) ~ — cos <pd{uj) 



nikiU^ z (0)u iyZ (z Q yJi(kix) 



(47) 



Two things can be noted. First, C vanishes anywhere above the source, (x = 0). The 
enhancement is destroyed because the dipolar nature of the source is in some sense "or- 
thogonal" to the detection of the displacement vector. Maximum enhancement actually 
occurs a fraction of a wavelength away from the source as shown in dashed line in Figure ||a. 
Secondly, the coherent enhancement around the source has a cos 2 structure, with </> the 
azimuthal angle between the dipole-axis of the source and the direction of detection. This 
"double- well" structure was observed by De Rosny, Tourin and Fink 31 ) . 

The coherent enhancement factor can be restored by a modification of the measurement. 
Suppose we measure the parameter d x u z (r,t). This measurement has the same symmetry 
as the dipolar source. We find for background and coherent enhancement, 

exp(— tot/Q) 1 



L(x,t) 



Dt 



~-d{u) 2 J^riikf \u i>z (0)\ J2rijk 2 \u jjZ 



C(x,t)^ eXp( -^ /g) x 



Dt 



d(u) 



J2nikiU i;Z (0)u iiZ (z o y 



J\{kjx) 



Jzfax) cos 2 



(48) 



For x = and z$ = we infer that L = C, i.e the maximal enhancement can now be reached. 
The plot of the restored enhancement factor as the function of the horizontal distance and 
as the function of the source depth are shown in solid line in the Figures [8|a, b. Note that 
the line profile is still not cylindrically symmetric, but depends on (j). 



D. Double-couple source at depth 

Seismic sources have successfully been modeled as two compensating couples (dipoles) 
TTJ . To facilitate observation of coherent backscattering with seismic waves we will here 



obtain the enhancement expected for such a source close to a free surface. In view of the 
complexity of the problem, we will restrict ourselves to a seismic plane that is oriented 
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parallel to the free surface where detection takes place. The depth of this plane is located 
at Zq. 

The force field of a double-couple source is described by a symmetric, off-diagonal seismic 
tensor. We assume that the two dipoles are orthogonal and along the axes x and y. The 
force field is then given by, 



f (r, u) = M{uo) {±d y + yd x ) S(r - r„) 



(49) 



with ro = (0,0,2;o). We can easily check that the mode representation of the source 
is S^k = u)M(u) [k x u ijy (z ) + kyU i>x (z )) We will assume that the measured parameter is 
d y 'U x i + d x iu y ', i.e. a certain horizontal component of the stress tensor; (x',y') are the 
coordinates in a frame that has been rotated over an angle j3 with (x, y) (see Figure |9|). The 
displacement vector of a mode (ik) can be expressed as, 



Uik(z) = {iii >z (z) z + uq (z) cos oti k + sin oci z x k | exp(ik 



(50) 



which introduces a new angle on independent on the direction k of propagation and on depth. 
Lamb waves have = whereas SH- waves have «j = |. We define as the angle between 
k and the x-axis, i.e. k = cos0x + sin0y. Finally, the angle \x orientates the direction of 
measurement x in the horizontal plane with respect to the source, (see Figure |9|). 
The incoherent background is calculated from, 



L ( X)t ) ^ exp(^/Q) M{uj) 2J2 f d 3 k [&u iky (0) + ^u ik , x ,(0)] 2 ImG ik l 



Dt 



10) 



X 



? / d3k ' + d y^'^ 2 Im G *' H 



(51) 



whereas the coherent enhancement follows from, 
C(M) ^ M(,f» 



J d 3 k [9 X /Uik, y '(0) + 9 y /u ik , x '(0)] [<9 x u ik , y (r ) + <9 y u ikiX (r )] ImG ik ( 



to, 



(52) 



These k-integrals can be evaluated straightforwardly and we simply quote the final result, 



(53) 



Here, wii denotes the complex amplitude of the horizontal component of the displacement 
vector. The coherent part is, 
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^2nikfui,\\(0)ui\\(z )* [cosf3J (kix) - cos q { J A (kix)} 



(54) 



with qi = 4/i + 3/9 for Lamb waves and ^ = 4/i + 3(3 + 7r for Sif waves. Since J\ term is very 
small, the line profile is almost isotropic around x = 0, independent on fi , and maximal 
for (3 — 0. The enhancement factor (L + C)/L is plotted in Figure pT0] a as a function of the 
horizontal distance for different source depth and in Figure |I~0fo as a function of the source 
depth z for a measurement on the top of the source. It is interesting to notice the relatively 
large second maximum of the coherent backscattering for a source at x — 0.7 X s away from 
the detector (x ~ 1.2km). 



VII. CONCLUSIONS AND OUTLOOK 



In this paper we have investigated multiple scattering of elastic waves for a model that 
is adapted to the needs of seismology of the Earth's crust: a two-dimensional solid plate 
with a thickness that is less than the mean free path of the waves. Contrary to other 
approaches, this model facilitates an exact treatment of the boundary condition at both 
sides of the plate, i.e. on the level of the elastic waves equation, and allows for fluctuations 
in the elastic constants that are depth-dependent. At the same time, we can describe the 
horizontal transport of waves, as well as the inter-mode mixing, by a generalized radiative 
transfer equation, that can be solved with conventional methods. Using this equation, 
we have investigated different aspects, such as surface detection, mode extinction times, 
equipartition, polarized sources at different depths in the plate and coherent backscattering. 
We believe that this study is an important step in the modelisation of seismic waves in the 
Earth's crust, but may also find applications in laboratory experiments. In future studies 
we will try to solve our equation numerically using Monte-Carlo methods. 

One final limitation of the present model has to be looked at in more detail. In this 
paper we have assumed a solid plate bounded by two ideal free surfaces without any leaks. 
The Earth's crust is much better described by one top free surface and a solid-solid interface 
at 30km in depth (the so-called Moho) overlying a high- velocity mantle. In previous studies 
we already suggested that the resulting energy leak into the mantle may be the origin of 
seismic coda |]23|j24]| . A future study should establish the relation between seismic coda and 
the individual quality factors of the modes. 
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TABLES 



N(u) 


3 


5 


13 


23 


43 


65 


85 


106 


D(u) 
D°°{iS) 


0.72 


0.56 


0.72 


0.77 


0.82 


0.84 


0.85 


0.85 



TABLE I. Ratio D(lv)/D°° (cj) as a function of the number of modes N(u) in the plate. 
D°°(uj) is the frequency-dependent diffusion constant for a 3D infinite medium, D(uj) is the fre- 
quency-dependent diffusion coefficient for our quasi-2D model with N(w) modes, with the same 
kind of disorder in A and fx. 
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FIGURES 

FIG. 1. Schematic plot of the dispersion law of the elastic Rayleigh-Lamb eigenmodes in a 
layer bounded by two free surfaces. Bold lines indicate symmetric branches, straight lines indi- 
cate anti-symmetric modes. Only modes of different symmetry are allowed to cross. The two 
dashed lines indicate the pure shear or pure compressional excitations. The surface Rayleigh waves 
propagate somewhat slower than pure S waves. 



FIG. 2. Extinction times for the different modes, calculated from Eq. (15), normalized by the 
mean free time of S waves in an infinite medium. The disorder is chosen to be uniform in the 
whole plate and the spatial correlation between the Lame coefficients is chosen equal. The plate 
thickness is H = 20.2A S , which has a number of N = 106 modes. 



FIG. 3. Eigenvalues of the matrix C defined in Eq. (2£) in the case of a plate thickness 
H = 20.2A S having a number of N = 106 modes normalized to the mean free time of 5-waves in 
an infinite medium. The disorder is uniform in the whole plate and the spatial correlation between 
all Lame coefficients is equal. 

FIG. 4. Energy distribution among the different modes for different sources. The total energy 
release is normalized to 1 and the plate thickness is H = 20.2A S with N = 106 modes. A) Isotropic 
explosion source at a depth \ s /3 from the free surface. B) Double couple source in the xy plan at 
a depth A s /3 from the free surface. C) Double couple source in the xz plan at a depth A s /3 from 
the free surface. The two Rayleigh modes are out of scale and carry half of the released energy. 
D) Double couple source in the xz plan at a depth 5A S from the free surface. 

FIG. 5. Prediction of the diffusion equation (^) for the time evolution of the energy for different 
sources. The time scale has been normalized to the mean free time of the S waves in an infinite 
medium. The plate thickness is H = 20.2A S , with N = 106 modes. A) and C) are predictions 
for the evolution of the energy for different modes, SH modes (Esh), Lamb modes without the 
Rayleigh modes (E( L _ R ^), and Rayleigh modes alone (Er) for respectively an isotropic explosion 
at a depth A s /3 from the free surface and a double couple source in the xz plan at depth 5A S from 
the free surface. B) and D) are predictions for the potential energy ratio E s /E p and the ratio 
H 2 /V 2 of the kinetic energies for the elastic displacement in different directions. 
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FIG. 6. Plot of the backscattering cone and enhancement factor for a monopolar source along 
the z axis. The normal component of the displacement field u z (0) is measured at the free surface 
and the plate thickness is H = 20.2A S which has N = 106 modes. A) Plot of the backscattering 
cone for different depths zq. B) Plot of the enhancement factor at x = as a function of the source 
depth zq. 

FIG. 7. Plot of the backscattering cone and enhancement factor for a isotropic explosion source. 
Both the divergence (solid line) and normal component of the field (dashed line) are measured. 
The plate thickness is H = 20.2A S which has N = 106 modes. A) Plot of the backscattering cone. 
B) Plot of the enhancement factor at x = as a function of the source depth zq. 

FIG. 8. Plot of the backscattering cone and enhancement factor for a dipolar source. Both 
d x u z (r,t) (solid line) and normal component u z (r,t) of the field (dashed line) are measured. The 
plate thickness is H = 20.2A S which has N = 106 modes. A) Plot of the backscattering cone. 
B) Plot of the enhancement factor at x = as a function of the source depth zq. 

FIG. 9. All angles involved in the measurement of the backscattering cone of a dislocation 
source at depth zq. See text for discussion, a« = for Lamb waves and ai = tt/2 for SH waves. 

FIG. 10. Plot of the backscattering cone and enhancement factor for a double-couple source 
with both its axes along the free surface. The orientation of the detection is such that (3 = and 
ix = 0. The plate thickness is H = 20.2A S which has N = 106 modes. A) Plot of the backscattering 
cone for different source depths zq = 0, zq = A s /3 and zq = X s /2. B) Plot of the enhancement 
factor at x = as a function of the source depth z$. 
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A) Monopolar Source 
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A) Dipolar Source 
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A) Double-Couple Source 
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